The wild silkworm Bombyx mandarina (Lepidoptera: Saturniidae) is considered one of the important pests whose larvae feed mainly on mulberry trees. Wild silkworm could be divided into two geographical groups based on the chromosomal characteristics, one is the Chinese Bombyx mandarina that mainly distributed in China and Russian Far East that number of chromosomes is the same as Bombyx mori (2n = 56), and the other is the Japanese Bombyx mandarina that mainly distributed in Japan and South Korea that number of chromosomes is different from Bombyx mori (2n = 54) (Arunkumar et al., 2006; Nakamura et al., 1999; Shimada et al., 1995). Research shown that silkworms may have been initially domesticated in China as tri-moulting lines, then subjected to independent spreads along the Silk Road that gave rise to the development of most local strains, and further improved for modern silk production in Japan and China, having descended from diverse ancestral sources (Lu et al., 2002; Pan et al., 2008; Sun et al., 2012; Xia et al., 2009; Xiang et al., 2018). Without artificial selection and live in the wild, wild silkworm preserved abundant genetic diversity, and shown significant differences in cocoon quality, disease resistance, stress tolerance, growth and development when compared with domestic silkworm, so it could be germplasm resources of genetic improvement of domestic silkworm (Fang et al., 2020; Fang et al., 2015; Xiang et al., 2013; Zhang et al., 2015). The origin and domestication of Bombyx mori and Bombyx mandarina are basic scientific issues in sericulture, and also involve in the inheritance of the silk culture. Therefore, clarifying the genetic relationship and geographical origination of silkworms based on the nuclear genome and mitogenome has important sense to sericulture development and scientific research (Goldsmith et al., 2005; Liu and Lu, 2018).

Mitochondria is a semi-autonomous organelle found in the majority of eukaryotic cells (Boore, 1999; Russell et al., 2020). Mitochondrial genomes (mtDNAs) are semi-independently replicating entities in mitochondria and also the only extra-nuclear genetic material in lots of animals. The mtDNA in most insects is typically a double-stranded circular molecule of 14~20 kb in length. The number of coding genes in mtDNA is thirty-seven that include twenty-two transfer RNA, two ribosome RNA and thirteen protein-coding genes (Boore, 1999; Pan et al., 2008). The mitogenome is widely used for population genetic structure, phylogenetic and phylogeography study for its characteristic of brief molecule structure, highly conservation, maternal inheritance, rapid evolutionary rate, less frequent reorganization and high copy number (Simon et al., 2006; Kuang et al., 2019; Chan et al., 2019; Min-Shan et al., 2018). The traditional methods to obtain mitogenome sequence are mainly based on sanger sequencing of PCR amplification. However, due to the lack of universal primer in many species, and the existence of special structures such as repetitive sequences and high AT content, it is difficult to obtain the full-length sequence of mtDNA by traditional methods. Therefore, high-throughput sequencing combined with mature bio-informatics methods have facilitate the assembly and annotation of mitogenome effectively, and become the most main technical way to obtaining the complete mitogenome (Kuang and Li, 2019).

Genome survey is an analysis to obtain bio-information such as heterozygote rate, GC content and genome size with low coverage high-throughput sequencing. It provides the basis for appropriate sequencing strategy and fine mapping of whole genome. For non-model organisms, the genome survey analysis is the basis for molecular mechanism study and genomic resources exploration. In this work, the genome size, GC content and heterozygote rate of Qin-ba wild silkworm were analyzed by Illumina NovaSeq high-throughput sequencing platform, and the mitogenome was assembled based on the sequencing data. The aim of this work is to provide a theoretical basis for whole genome sequencing, and the cluster analysis of silkworms from different geographical origins could contributes to a better comprehension of the genetic diversity and domestication of silkworm. This work could offer theoretical references for utilization of wild silkworm and inheritance of silk culture.

1 Results and Conclusion

1.1 Sequencing data statistics

A total of 26.8 Gb sequence data were generated from the small-insert (300 bp) library of Bombyx mandarina. A total of 25.8 Gb clean bases were generated after the sequence data was filtered, low-quality and duplication sequences were trimmed, with 92% Q20 bases and 38.0% GC content which was used for assembly.

1.2 prediction of genome size, heterozygosity ratio

The estimation of genome size and heterozygosity ratio were based on the K-mer analysis (Figure 1). Statistical analysis showed that the total number of K-mers was 22 366 412 989. The results shown that the 17-mer frequency distribution curve exhibited two peaks at depths of 24 and 49, respectively. Using the formula of genome size = total K-mer number / peak depth, the genome size of Qin-ba wild silkworm was estimated to be 456.5 MB. The first peak observed at 1/2 of peak depth displayed a high level of heterozygosity for this wild silkworm sample. Simulation analysis using the Arabidopsis thaliana genome revealed that it had a 1.94% heterozygosity rate, represent that the genome belonging to the complex genome with higher heterozygosis rate.

Figure 1 17-mer distribution and heterozygosity simulation curve of Bombyx mandarina

1.3 Preliminary genome assembly and GC content-sequencing depth analysis

SOAP denovo program was employed to preliminary genome assembly after filtered and removed low quality data reads (Table 1). The results shown that the software produced a contig with the N50 of 587 bp, the longest contig length of 23.641 kb, and the total length of 540.745 MB, while sequence with a scaffold N50 length of 1792 bp, the longest scaffold length of 64.381 kb, and the total length of 649.282 Mb.

Table 1 Statistics of assembly sequences

The statistics of GC content and read depth were shown as Figure 2. It represents that most plots scatted in the area of GC 30~40% and depth 40~60, implied that there is no contaminant in the sequencing reads, and the heterozygosity result the read-depth clusters divided into two obvious layers. Ten thousand scaffolds were random selected to align to the NCBI Nucleotide (Nt) database. Its shown that 51.17% of them matched with Bombyx mori, while 47.14% of them acted no-hits. The results of blast showed that these samples do not contain exogenous contamination.

Figure 2 Statistics of GC content and read depth

1.4 Assembly and annotation of Mitochondrial genome

The full-length mitogenome of Bombyx mandarina was assembled from the high-throughput sequencing data, shown as Figure 3. Its shown that the mitogenome size is 15 662 bp with 37 genes in it. The number of genes for protein coding, transfer RNA and ribosomal RNA are 13, 22 and 2 respectively, while the length of transfer RNA genes variated from 61~71 bp. The protein coding genes include one cytochrome b apoenzyme (CYTB), two ATP synthase subunits (ATP6 and ATP8), three cytochrome oxidase subunits (COX1, COX2 and COX3) and seven NADH dehydrogenase subunits (ND1, ND2, ND3, ND4L, ND4, ND5 and ND6). Twenty-five intergenic- and five overlapping- regions distributed in the mitogenome of Qin-ba Bombyx mandarina (Table 2). The gene order of wild silkworm was the same as that in fruit fly and domestic silkworm, and no re-arrangement observed.

Figure 3 Structure of the mitochondrial genome of Bombyx mandarina

Table 2 Structure of the mitochondrial genome of Bombyx mandarina

1.5 Phylogenetic analysis

For cluster analysis of wild- and domestic- silkworm, the mitogenomes of silkworms were downloaded from GenBank to constructed the phylogenetic tree, and Drosophila melanogaster, Antheraea pernyi and Rondotia menciana were set as the out group (Figure 4). The result shown that the wild silkworms could be divided into Japanese, southern or northern of china sub groups according to their origin geographical areas. Mitogenome of wild silkworm that from Shandong, Liaoning and Shaanxi belong to the sub group of northern china. The wild silkworms from northern china have the closer relationship with domestic silkworm than that of Japanese and southern china.

Figure 4 Phylogenetic tree of Bombyx mori and Bombyx mandarina based on the mitochondrial genome. The “Bman” represent the Bombyx mandarina, “Bm” represent the Bombyx mori. bootstrap values were marked on the branches

2 Discussion

For non-model organism without genomic data, study on the genomic features could pave the foundation for molecular mechanism research and gene resource utilization. Flow cytometry and Feulgen spectrophotometry are the common approaches for genomic size estimation, while high throughput sequencing is another faster route to complete genomic survey. In this work, we carried out genome survey of wild silkworm from Ankang city, Shaanxi province, via the high-throughput sequencing. The results shown that the estimated genome size is 456.5 Mb, and the heterozygosity ratio is 1.94%. The preliminary genome assembly produced a contig with the N50 of 587 bp, the longest contig length of 23.641 kb, and the total length of 540.745 MB, while scaffold N50 length of 1792 bp, the longest scaffold length of 64.381 kb, and the total length of 649.282 Mb. The complete mitogenome length is 15662 bp with 37 genes in it, no gene arrangement observed. As phylogenetic analysis based on the mitogenomes, the china wild silkworm could be divided into northern and southern geographic sub groups, and the northern sub group has the closest relationship with domestic silkworm.

In general, genome with heterozygosity ratio up than 0.8% and repeat sequence ration outweigh 60% could be considered as complex genome, then it is hardly to obtain high quality genome through the conventional sequencing and assembly technology (Gao et al., 2018). Our work represents that genome of Qin-ba wild silkworm belongs to complex genome, and the results of preliminary genome assembly through Illumina sequencing is not ideal. Studies shown that it is hardly for the second- or the third-generation sequencing to overcome the repeats in genomes, while high-throughput/resolution chromosome conformation capture technology (Hi-C) can depict three-dimensional global chromatin interactions across eukaryotic genome. Therefore, integrating the third-generation sequencing technology with the Hi-C could be great helpful to obtain high quality genome of wild silkworm.

The wild silkworm Bombyx mandarina is the wild ancestor of the domestic silkworm Bombyx mori, but there remain single- or multi- geographical origin hypothesis of domestication. Chen DB revealed that Chinese B. mandarina populations represented two genetically distinctive subtypes in line with the geographic boundary of northern and southern China based on the mitogenome analysis, and the true wild ancestor of domestic silkworm is northern Chinese Bombyx mandarina, rather than southern Chinese Bombyx mandarina (Chen et al., 2019). Our work supported this notion, and implied that the Qin-ba mountain area could be one geography region of the origin of domestic silkworm.

With of NGS technology emergency and decreased sequencing cost, functional genomics research has been promoted greatly as more and more insect genomes have been analyzed (Zhao, 2016). The genomic research of wild silkworm could be conductive to the development of pest control. Otherwise, the comparative genomics study in wild- and domestic- silkworm could result in new discovery in domestic mechanisms and gene resources, and facilitate the new germplasm creation and wild resource utilization. In this work, we survived the genome of Qin-ba wild silkworm by using high-throughput sequencing technology, and analyzed phylogenetic relationships of silkworm based on the mitogenomes, could provide a scientific basis for the future research on the whole genome mapping of wild silkworm.

3 Materials and Methods

3.1 Experimental materials

The wild silkworm was collected from Raofeng town, Shiquan county, Ankang City, and raising in Shaanxi Key Laboratory of Sericulture more than six years as conventional mulberry culture. One female individual was used to extracted genome, and then conserved under ultralow temperature for use.

3.2 Sample extraction and Sequencing data generation

Genomic DNA was extracted from the pupal sample by using DNeasy Blood & Tissue Kit (Qiagen, Germany) following the manufacturer’s protocol. UV spectrophotometry was used to measure the concentration and purity the DNA. Agarose gel electrophoresis was used to determine the integrity of the template, and then sheared randomly into small fragments by ultrasonic wave. Two paired-end libraries with an insert size of 300 base pairs (bp) were constructed from fragmented random genomic DNA following the Illumina manufacturer’s instructions and sequenced by Illumina NovaSeq sequencing platform. The clean reads were obtained by removing reads containing adapter and low-quality reads from raw data with NGS-QC-Generator Toolkit. All the downstream analyses were based on these clean data.

3.3 Prediction of genome size, heterozygosity ratio and Preliminary genome assembly

K-mer analysis was used to predict the genome size and heterozygosity ratio before genome assembly. A K-value of 17 was used for the prediction, analysis, and iterative selection of 17-bp base sequences from the clean reads. SOAPdenovo software was used to carry out preliminary genome assembly after the short tips and low-quality sequences were filtered.

3.4 Analysis of GC content and GC-depth statistics

The reads obtained from sequencing all the libraries were aligned to the initially obtained contigs. The filtered reads were aligned to this assembled sequence using SOAP to obtain the base depth. A window size of 10 kb was used for non-repetitive advancement in the sequence and calculation of the mean depth and GC content of every window to generate a GC depth plot to examine whether there was significant GC bias or bacterial contamination in the sequences.

3.5 Assembly and annotation of mitochondrial genome

The mitogenome was built using the NOVOPlasty pipeline v. 2.6.3 with default setting, and the sequence with GenBank no. MG604734.1 set as the reference sequence. The newly assembled mitochondrial genome was annotated and visualized in the GeSeq web server using the invertebrate genetic code. The mitogenome of this work were uploaded to GenBank with the accession no. MN400656.

3.6 Phylogenetic analysis

The mitogenomes of 25 domestic silkworms and 18 wild silkworms were downloaded from GenBank. The phylogenetic tree based on neighbor-joining algorithm, and genetic distance analysis based on Kimura 2-parameter model were conducted by MEGA-X software. The mitogenomes of Drosophila melanogaster, Antheraea pernyi and Rondotia menciana were set as the out group.

Authors’ contributions

MG was project leader as performed the experimental design, data analysis and article wrote. WRX participated in article wrote and experimental design. CQ participated in experiment operation and data analysis. All authors read and approved the final manuscript.

Acknowledgments

This work was supported by the Key R & D plan of Shaanxi Province (No.2020NY-014); Opening Foundation of State Key Laboratory of Silkworm Genome Biology (No. SKLSQB1819-4); Educational Commission of Shaanxi Province of China (No.18JS001 and No.18JS003).

References

Arunkumar K.P., Metta M., and Nagaraju J., 2006, Molecular phylogeny of silkmoths reveals the origin of domesticated silkmoth, Bombyx mori from Chinese Bombyx mandarina and paternal inheritance of Antheraea proylei mitochondrial DNA, Molecular Phylogenetics and Evolution, 40(2): 419-427

https://doi.org/10.1016/j.ympev.2006.02.023

PMid:16644243

Boore J.L., 1999, Animal mitochondrial genomes, Nucleic Acids Research, 27(8): 1767-1780

https://doi.org/10.1093/nar/27.8.1767

PMid:10101183 PMCid:PMC148383

Chan E.K.F., Timmermann A., Baldi B.F., Moore A.E., Lyons R.J., Lee S.S., Kalsbeek A.M.F., Petersen D.C., Rautenbach H.,Fortsch H.E.A., Bornman M.S.R., and Hayes V.M., 2019, Human origins in a southern African palaeo-wetland and first migrations, Nature, 575(7781): 185-189

https://doi.org/10.1038/s41586-019-1714-1

PMid:31659339

Chen D.B., Zhang R.S., Bian H.X., Li Q., Xia R.X., Li Y.P., Liu Y.Q., and Lu C., 2019, Comparative mitochondrial genomes provide new insights into the true wild progenitor and origin of domestic silkworm Bombyx mori, International Journal of Biological Macromolecules, 131: 176-183

https://doi.org/10.1016/j.ijbiomac.2019.03.002

PMid:30836184

Fang S.M., Hu B.L., Zhou Q.Z., Yu Q.Y., and Zhang Z., 2015, Comparative analysis of the silk gland transcriptomes between the domestic and wild silkworms, BMC Genomics, 16(1): 60-72

https://doi.org/10.1186/s12864-015-1287-9

PMid:25887670 PMCid:PMC4328555

Fang S.M., Zhou Q.Z., Yu Q.Y., and Zhang Z., 2020, Genetic and genomic analysis for cocoon yield traits in silkworm, Scientific Reports, 10(1): 5682-5693

https://doi.org/10.1038/s41598-020-62507-9

PMid:32231221 PMCid:PMC7105477

Gao S.H., Yu H.Y., Wu S.Y., Wang S., Geng J.N., Luo Y.F., and Hu S.N., 2018, Advances of sequencing and assembling technologies for complex genomes, Hereditas, 40(11): 944-963

Goldsmith M.R., Shimada T., and Abe H., 2005, The genetics and genomics of the silkworm, Bombyx mori, Annual Review of Entomology, 50: 71-100

https://doi.org/10.1146/annurev.ento.50.071803.130456

PMid:15355234

Kuang W.M., and Li Y., 2019, Mitogenome assembly strategies and software applications in the genome era, Hereditas, 41(11): 979-993

Kuang W.M., Ming C., Li H., Wu H., Frantz L., Roos C., Zhang Y., Zhang C., Jia T., Yang J., and Yu L., 2019, The Origin and Population History of the Endangered Golden Snub-Nosed Monkey (Rhinopithecus roxellana), Molecular Biology and Evolution, 36: 487-499

https://doi.org/10.1093/molbev/msy220

PMid:30481341

Kumar S., Stecher G., Li M., Knyaz C., and Tamura K., 2018, MEGA X: Molecular Evolutionary Genetics Analysis across computing platforms, Molecular Biology & Evolution, 35(6): 1547-1549

https://doi.org/10.1093/molbev/msy096

PMid:29722887 PMCid:PMC5967553

Liu Y.Q., and Lu C., 2018, Research advances in origin and evolution of domestic Silkworm, Sericulture Science, 44(3): 353-358

Lu C., Yu H.S., and Xiang Z.H., 2002, Molecular Systematic Studies on Chinese Mandarina Silkworm (Bombyx mandarina M. ) and Domestic Silkworm (Bombyx mori L.), Agricultural Sciences in China, 1(3): 349-358

Min-Shan K.A., Zhang Y., Yang M.A., Hu Y., Cao P., Feng X., Zhang L., Wei F., and Fu Q., 2018, Mitochondrial genome of a 22,000-year-old giant panda from southern China reveals a new panda lineage, Current Biology, 28(12): R693-R694

https://doi.org/10.1016/j.cub.2018.05.008

PMid:29920259

Nakamura T., Banno Y., Nakada T., Nho S.K., Xu M.K., Ueda K., Kawarabata T., Kawaguchi Y., and Koga K., 1999, Geographic dimorphism of the wild silkworm, Bombyx mandarina, in the chromosome number and the occurrence of a retroposon-like insertion in the arylphorin gene, Genome, 42(6): 1117-1120

https://doi.org/10.1139/gen-42-6-1117

PMid:10659778

Nicolas D., Patrick M., and Guillaume S., 2016, NOVOPlasty: de novo assembly of organelle genomes from whole genome data, Nucleic Acids Research, 45(4), e18, doi: 10.1093/nar/gkw955

https://doi.org/10.1093/nar/gkw955

PMid:28204566 PMCid:PMC5389512

Pan M.H., Yu Q.Y. Xia Y.L., Dai F.Y., Liu Y.Q., Lu C., Zhang Z., and Xiang Z.H., 2008, The Characteristics of mitochondrial genome of Chinese wild silkworm (Bombyx mandarina), Science China, 36(8): 751-759

https://doi.org/10.1007/s11427-008-0097-6

PMid:18677597

Russell O.M., Gorman G.S., Lightowlers R.N., and Turnbull D.M., 2020, Mitochondrial Diseases: Hope for the Future, Cell, 181(1): 168-188

https://doi.org/10.1016/j.cell.2020.02.051

PMid:32220313

Shimada T., Kurimoto Y., and Kobayashi M., 1995, Phylogenetic Relationship of Silkmoths Inferred from Sequence Data of the Arylphorin Gene, Molecular Phylogenetics & Evolution, 4(3): 223-234

https://doi.org/10.1006/mpev.1995.1021

PMid:8845960

Simon C., Buckley T.R., Frati F., Stewart J.B., and Beckenbach A.T., 2006, Incorporating Molecular Evolution into Phylogenetic Analysis, and a New Compilation of Conserved Polymerase Chain Reaction Primers for Animal Mitochondrial DNA, Annual Review of Ecology, Evolution, and Systematic, 37: 545-579

https://doi.org/10.1146/annurev.ecolsys.37.091305.110018

Sun W., Yu H.S., Shen Y.H., Banno Y., Xiang Z.H., and Zhang Z., 2012, Phylogeny and evolutionary history of the silkworm, Science China Life Sciences, 55(6): 483-496

https://doi.org/10.1007/s11427-012-4334-7

PMid:22744178

 Tillich M., Lehwark P., Pellizzer T., Ulbricht-Jones E.S., Fischer A., Bock R., and Greiner S., 2017, GeSeq - versatile and accurate annotation of organelle genomes, Nucleic Acids Research, 45(w1): W6-W11

https://doi.org/10.1093/nar/gkx391

PMid:28486635 PMCid:PMC5570176

Xiang H., Li X., Dai F., Xu X., Tan A., Chen L., Zhang G., Ding Y., Li Q., Lian J., Willden A., Guo Q., Xia Q., Wang J., and Wang W., 2013, Comparative methylomics between domesticated and wild silkworms implies possible epigenetic influences on silkworm domestication, BMC Genomics, 14: 646-656

https://doi.org/10.1186/1471-2164-14-646

PMid:24059350 PMCid:PMC3852238

Xia Q., Guo Y., Zhang Z., Li D., Xuan Z., Li Z., Dai F., Li Y., Cheng D., Li R., Cheng T., Jiang T., Becquet C., Xu X., Liu C., Zha X., Fan W., Lin Y., Shen Y., Jiang L., Jensen J., Hellmann I., Tang S., Zhao P., Xu H., Yu C., Zhang G., Li J., Cao J., Liu S., He N., Zhou Y., Liu H., Zhao J., Ye C., Du Z., Pan G., Zhao A., Shao H., Zeng W., Wu P., Li C., Pan M., Li J., Yin X., Li D., Wang J., Zheng H., Wang W., Zhang X., Li S., Yang H., Lu C., Nielsen R., Zhou Z., Wang J., Xiang Z., and Wang J., 2009, Complete resequencing of 40 genomes reveals domestication events and genes in silkworm (Bombyx), Science, 326(5951): 433-436

https://doi.org/10.1126/science.1176620

PMid:19713493 PMCid:PMC3951477

Xiang H., Liu X., Li M., Zhu Y., Wang L., Cui Y., Liu L., Fang G., Qian H., Xu A., Wang W., and Zhan S., 2018, The evolutionary road from wild moth to domestic silkworm, Nature Ecology & Evolution, 2(8): 1268-1279

https://doi.org/10.1038/s41559-018-0593-4

PMid:29967484

Zhang X.T., Nie M.Y., Zhao Q., Wu Y.Q., Wang G.H, and Xia Q.Y, 2015, Genome-wide patterns of genetic variation among silkworms, Molecular Genetics and Genomics, 290(4): 1575-1587

https://doi.org/10.1007/s00438-015-1017-7

PMid:25749967

Zhao X.F., 2016, Basic scientific questions and molecular biology techniques in entomology, Acta Entomologica Sinica, 59(8): 896-905